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^n^j , We present a review of our recent work in extending the successful dy- 

namical mean-field theory from the equilibrium case to nonequilibrium 
' cases. In particular, we focus on the problem of turning on a spatially 

. uniform, but possibly time varying, electric field (neglecting all magnetic 

field effects). We show how to work with a manifestly gauge-invariant 
formalism, and compare numerical calculations from a transient-response 
formalism to different types of approximate treatments, including the 
semiclassical Boltzmann equation and perturbation theory in the inter- 
action. In this review, we solve the nonequilibrium problem for the 
Falicov-Kimball model, which is the simplest many-body model and the 
easiest problem to illustrate the nonequilibrium behavior in both dif- 
fusive metals and Mott insulators. Due to space restrictions, we as- 
sume the reader already has some familiarity both with the Kadanoff- 
Baym-Keldysh nonequilibrium formalism and with equilibrium dynam- 
ical mean-field theory; we provide a guide to the literature where addi- 



X 

H 

' tional details can be found. 



1.1. Introduction 

Dynamical mean-field theory (DMFT) was introduced in 1989 1 shortly af- 
ter Metzner and Vollhardt 2 proposed scaling the hopping matrix clement 
as the inverse square root of the spatial dimension to achieve a nontriv- 
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ial limit where the many-body dynamics are local. Since then, the field 
has blossomed to the point where nearly all model many-body problems 
have now been solved, 3 and much recent work has focused on applying 
DMFT principles to real materials calculations. 4 Little work has empha- 
sized nonequilibrium aspects of the many-body problem, where the strongly 
correlated system is driven by an external field that can possibly sustain 
a nonequilibrium steady state. In this contribution, we will review recent 
work that has been completed on expanding DMFT approaches into the 
nonequilibrium realm. Wc will show how to work with so-called gauge- 
invariant Green functions 5 to illustrate that one can carry out calculations 
in a form that manifestly is independent of the gauge chosen to describe 
the driving fields. This approach is different from our previously published 
work, where we worked solely with Green functions in the Hamiltonian 
gauge (where the scalar potential vanishes). 

We examine the problem of strongly correlated electrons driven by a 
spatially uniform electric field in the limit of infinite dimensions, 6-8 where 
DMFT can be applied to solve the problem exactly. In infinite dimen- 
sions, the self-energy of the electrons is local, and the lattice problem can 
be mapped onto the problem of an impurity coupled to an effective time- 
dependent field (which is adjusted so that the impurity Green function and 
the local Green function on the lattice are identical). The impurity prob- 
lem in the dynamical mean field can be solved exactly for many different 
cases. In equilibrium, a large number of strongly correlated models have 
been solved in infinite dimensions, like the Falicov-Kimball model, 1,10 ' 11 the 
Hubbard model, 12-14 the periodic Anderson model, 15 ' 16 and the Holstein 
model, 17 ' 18 (for a reviews, see Refs. 3 and 19). Recently, there has been a 
significant effort in combining DMFT with density functional theory (DFT) 
to describe properties of real materials when DFT is insufficient to properly 
describe the electron-electron interactions (see Ref. 4 for a review). It is 
now generally believed that DMFT is a good approximation to the many- 
body problem in three dimensions, and it can accurately describe strong 
electron-electron correlation effects in bulk systems. 

The first attempt to employ DMFT to describe nonequilibrium prop- 
erties of a strongly correlated model was made by Schmidt and Monien in 
Ref. 20, where they studied the spectral properties of the Hubbard model in 
the presence of a timc-dependent chemical potential by using iterated per- 
turbation theory (PT). Recently, we have developed a generalized nonequi- 
librium DMFT formalism to study the response of correlated electrons to a 
spatially uniform time-dependent electric field and applied that formalism 
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to the Falicov-Kimball model. 6-8 The Falicov-Kimball model, 21 is the sim- 
plest model for strongly correlated electrons that demonstrates long range 
order and undergoes a metal-to-Mott-insulator transition. It consists of two 
kinds of electrons: conducting c-electrons and localized /-electrons, which 
interact through an on-site Coulomb repulsion. The model was introduced 
to describe valence- change and metal-insulator transitions 21 in rare-earth 
and transition-metal compounds. It was reinvented as a model to describe 
crystal formation 22 resulting from the Pauli exclusion principle. DMFT 
was actually developed with the original solution of the Falicov-Kimball 
model in infinite dimensions 1,10,11 and now there is an almost complete 
understanding of its general properties (for a review, see Ref. 19). We 
extended the equilibrium formalism to the nonequilibrium case, where we 
numerically solved a system of the equations for the Green function and 
self-energy defined on a complex time contour (see Fig. II. ip by using the 
Kadanoff-Baym-Keldysh nonequilibrium Green function formalism. 23,24 

In this review, we summarize the successes of recent work to generalize 
DMFT to nonequilibrium problems with a focus on solutions of the spinless 
Falicov-Kimball model on an infinite-dimensional hypercubic lattice in the 
presence of an external time-dependent electric field. There are many in- 
teresting and surprising results which differ from semiclassical predictions 
(such as those made from the Boltzmann equation solution). In addition 
to the exact solutions, we also present results for the noninteracting case 
and for the case of second-order perturbation theory in the interaction. In 
particular, we analyze the limitations of the perturbation theory approxi- 
mation, especially in studying (long-time) steady-state behavior. 

1.2. General nonequilibrium formalism 

The nonequilibrium properties of a quantum many-particle system can be 
studied by calculating the contour-ordered Green function in momentum 
space: 



defined on the complex time-contour presented in Fig. 11.11 (see, for ex- 
ample, Ref. 25). Since the system is initially in equilibrium, the ther- 
mal average in Eq. (|1.1[) is performed with the equilibrium density matrix 
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Fig. 1.1. The complex Kadanoff-Baym-Kcldysh time contour for the two-time Green 
functions in the nonequilibrium case. The time increases from —t ma x (left point on 
the top branch) along the contour to t max then decreases back to —t max and then runs 
parallel to the imaginary axis to —t ma x — if). We consider the situation when the electric 
field is turned on at t = 0, so the vector potential is nonzero for t > 0. We assume that 
both t\ and t2 lie somewhere on the contour. 

exp[— f3H(— t max )]/Trexp[— (3H{— t max )] with respect to the initial Hamil- 
tonian H(—t max ) with vanishing electric field (the symbol (3 = 1/T is the 
inverse temperature). The operator indices H and / in Eq. (|1.1|) stand 
for the Heisenberg and Interaction representations, respectively. In this 
formalism, familiar quantum many-body techniques derived in equilibrium 
can also be used in the nonequilibrium case, except that now the time or- 
dering T c of the operators is along the complex contour. In particular, 
the Schwinger-Dyson equation, which connects the contour-ordered Green 
function with the electron self-energy ££.(£i,i2), remains valid: 



where the matrix product of the continuous matrix operators is accom- 
plished by line integrals over the contour. 

In DMFT, we work with the local Green function, which is found by 
summing the momentum-dependent Green function over all momenta. We 
then map the many-body problem on the lattice to an impurity problem, 
but in a dynamical mean field that mimics the hopping of electrons onto and 
off of the given site. It turns out that one needs the full freedom available 
with the three-branch contour to find the proper dynamical mean field to 
map the impurity onto the lattice. Hence, our approach will work with 
the less common Green functions on the three-branch contour, as opposed 
to a simpler two-branch contour, which we work with when we discuss 
the perturbative approach on the lattice. One can find the time-ordered, 
anti-time-ordered, lesser, greater, retarded, advanced and thermal Green 
functions on this contour. 26 For example, the retarded Green function, 




(1.2) 
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which is related to the density of quantum states, is 

GQ(h,t2) = -ie(h - h){{c^ H {t 1 ),cl H {t 2 )} + ), (1.3) 

where the braces indicate the anticommutator of the two operators, and the 
lesser Green function, which is related to how the electrons are distributed 
amongst the quantum states, satisfies 

G^(ii,ta)=*(cL r (t2)c kH (*i)>. (1-4) 
Both of these functions can be extracted from G£. 

1.3. Nonequilibrium dynamical mean-field theory for the 
Falicov-Kimball model 

The spinless Falicov-Kimball model 21 consists of two kinds of electrons: 
conduction c-electrons and localized /-electrons. They interact with each 
other through an on-site Coulomb repulsion U . The model Hamiltonian 
has the following form in the absence of any external fields: 

n = - e uAtj - mE c \°i - n E fifi + u E fifi&i, (i-s) 

(ij) i i i 

where tij = t is the nearest-neighbor hopping matrix clement for the c- 
electrons, \x and /j,f are the chemical potentials of c- and /-electrons, corre- 
spondingly. Due to the Pauli principle, there is no local cc- and //-electron 
interaction in the spinless case. The Hamiltonian in Eq. (| 1 . 5j) can also be 
regarded as an approximation to the spin s = 1/2 Hubbard model, where 
spin-up (o) electrons move in a frozen background of the localized spin- 
down (/-electrons). We consider the problem on the infinite-dimensional 
(d — * oo ) hypercubic lattice at half-filling, when the particle densities of 
the c- and /-electrons are equal to 0.5. In this limit, the hopping pa- 
rameter is renormalizcd in the following way: 2 t = t* j1\fd. In the limit 
of infinite dimensions, one can solve the equilibrium problem for the con- 
duction electrons exactly at any temperature, particle concentration and 
Coulomb repulsion. The key simplification, which allows one to obtain 
the exact solution as d — > oo, comes from the fact that the electron self- 
energy is momentum-independent. 1,27 Although that original work was 
performed in equilibrium, Langreth's rules 28 guarantee that it also holds 
for the nonequilibrium case. 

Nowadays, most of the equilibrium properties of the model, including 
the phase diagram, are well known (see Ref. 19). In particular, the model 
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demonstrates a Mott transition when n c + rif = 1 at some critical value 
of the Coulomb repulsion, 29 which depends on the particular value of n c 
(rif is equal to 1 — n c in this case). In the insulating phase, the density 
of states A(cu) is not equal to zero for frequencies inside the "gap region" , 
but is exponentially suppressed, except for uo = 0. Therefore, the density of 
states actually demonstrates a pseudogap in the insulating phase, which is 
an artifact of the fact that the infinite-dimensional hypercubic lattice has 
a Gaussian density of states for the noninteracting problem, which does 
not have a finite bandwidth. Another important feature is the behavior 
of the imaginary part of the self-energy for frequencies close to zero in the 
"metallic" phase: ImE(w) ~ — c + e'eu 2 (c and d > and independent 
of temperature), which differs from the standard Fermi liquid behavior 

ImE(cj) a(T) - bu; 2 (a and b > and a(T) -► as T -► 0). This means 

that there are no long-lived Fermi liquid quasiparticles in the model. 

We are interested in the case when the system is coupled to an exter- 
nal electromagnetic field E(r,t). This field can be expressed by a scalar 
potential <p(r, t) and by a vector potential A(r, t) in the following way: 

E(r,t) = -V^(r, f )-i^^. (1.6) 

We assume that the electric field is spatially uniform and choose the tem- 
poral or Hamiltonian gauge for the electric field: ip(r, t) = 0. In this case, 
the electric field is introduced into the Hamiltonian by means of the Peierls 
substitution for the hopping matrix: 30 



tij ^ tij exp 



le 
he 



/ A(r,t)dr 



tij exp 



A(i) • (Ri - Rj) 



(1.7) 

where the last formula holds for a spatially uniform field where we take 
A(i) = — Eci#(i) for a uniform field turned on at t = 0. 

We assume that it is safe to neglect magnetic field effects, because the 
electric field varies slow enough in time (recall Maxwell's equations say that 
a time- varying electric field creates a time varying magnetic field). This 
approximation is valid when the electric field is smooth enough in time 
that the magnetic fields can be ignored. Another way of describing this is 
that we assume our electric field is always spatially uniform, even though 
it has a time dependence, which is not precisely a solution of Maxwell's 
equations, but is approximately so. 

The electric field introduced into the Hamiltonian Eq. (|1.5p results in 
a time-dependent shift of the momentum in the free electron dispersion 
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eA(t) 

he 



-It cos 



1=1 



a k 



he 



(1.8) 



It is convenient to consider the case, when the electric field lies along the 
elementary cell diagonal: 31 



A(i)=A(t)(l,l,... J l). 
In this case, the free electron spectrum 



(1.9) 



cos (2^2) «(k) + si, , (k), (1.10) 



depends on only two energy functions: 

e(k) = -2t^2cos(ak l ) 



and 



e(k) = -2i^sin(a£;'). 



(1.11) 



(1.12) 



Of course, when the field vanishes, the energy spectra in Eq. p.!0[) reduces 
to the standard spectra in Eq. (|1.11[) for free electrons on the hypercu- 
bic lattice. In the limit of an infinite dimensional hypercubic lattice, one 
can calculate the joint density of states for the two energy functions in 
Eqs. ([01]) and (fLT2"j) . 32 



P2(e,e) 



1 



7rf*V 



■ exp 



e 



s 

^2 



(1.13) 



Below we use atomic units, putting all fundamental constants, except the 
electron charge e, to be equal to one: a = H = c = t* = l. 

To solve the problem of the response of the conduction electrons to an 
external electric field, we use a generalized nonequilibrium DMFT formal- 
ism. 8 The electron Green functions and self-energies are functions of two 
time arguments defined on the complex time-contour in Fig. 11.11 Since the 
action for the Falicov-Kimball model is quadratic in the conduction elec- 
trons, the Feynman path integral over the Kadanoff-Baym-Keldysh contour 
can be expressed by the determinant of a continuous matrix operator with 
arguments defined on the contour. Because the concentration of localized 
particles on each site is conserved, one can calculate the trace over the 
fcrmionic variables. It is possible to show that the self-energy remains local 
in the limit of infinite dimensions in the presence of a field; start with the 
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equilibrium perturbation theory expansion for the self-energy 2 and then 
apply Langreth's rules 28 to the self-energy diagrams, which say that every 
nonequilibrium diagram is obtained from a corresponding equilibrium dia- 
gram, with the time variables now defined on the Kadanoff-Baym-Keldysh 
contour. 

The generalized system of nonequilibrium DMFT equations for the con- 
tour ordered Green function G(£i,£ 2 ), self-energy £(£i,£2) and an effective 
dynamical mean-field A(£i,£2) can be written in analogy with the equilib- 
rium case 12 as follows: 



G(£ 1; £ 2 ) = (1 - u/i)G„(ti,t2) + «Ji[G&L,(/i - U) - A]- 1 (« 1) t 2 )(1.17) 



where G^ (£1, £2) is the noninteracting electron Green function in the pres- 
ence of an external time-dependent electric field, which can be calculated 
analytically (see below) and Goimp(ti, £2; is the free impurity Green func- 
tion in a chemical potential fi. The symbol w\ is the average number of the 
/-electrons per site. In our case, Wi = 1/2. 

The momentum summation in Eq. (|1.14| can be performed by intro- 
ducing the two energy functions Eqs. fll.lljl and (|1.12|) and using the joint 
density of states in Eq. (|1.13|) : J2k ^k = / de J dep2(e, e)F^ tS whenever the 
summand depends on momentum only through the two energy functions. 
The system of equations (|1.14j) - (|1.17[) formally resembles the corresponding 
system in the equilibrium case, except now we have to work with a two- 
time formalism on the contour, rather than being able to Fourier transform 
the relative time to a frequency. And, because we are working with the 
contour-ordered Green functions, which depend on the distribution of elec- 
trons, we need to be careful to treat how the chemical potential is shifted 
by U when we perform the trace over the /-electrons. 

The system of equations (I1.14[) - (I1.17I) can be solved by iteration as fol- 
lows. One starts with an initial self-energy matrix, for example the equi- 
librium self-energy. Substitution of this function into Eq. (|1.14[) gives the 
Green function. Then, from Eq. (|1.16[) one can find the effective dynamical 
mean-field A(ii,i 2 ), which allows one to find the new value for the Green 
function G(t±, £2) from Eq. (|1.17[) . After that, one finds the new self-energy 
£(£i,£ 2 ) from the impurity Dyson equation and the dynamical mean field. 




(1.14) 



k 



Go(tx, £2) 
A(£i,£ 2 ) 



^0imp(*l'*2;M) _ G 1 (£i,£ 2 ), 



(1.15) 
(1.16) 
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The calculations are repeated until the difference between the old and new 
values for the self-energy £(ti,t2) are smaller than some desired precision 
(usually 10 -6 in relative error). 

In practice, to solve this system numerically, one needs to discretize the 
complex time contour Fig. (jl.ljl with some step At along the real axis and 
At along the imaginary axis. In this case, the functions in Eqs. (II. 14[) - 
(I1.17P become general complex square matrices. In order to study the long 
time behavior, one needs to choose the value of t max large enough. The 
precision of the solution strongly depends on the value of the discretization 
step At, which must be small enough. Therefore, in order to get a precise 
long time solution it is necessary to use large complex square matrices 
in Eqs. (|1.14|) - (|1.17|) . This causes some constraints connected with the 
machine memory and the computational time. In our calculations, we used 
the time step At ranging from 0.1 to 0.0167 and matrices up to order 
4900 x 4900. The two-energy integration in Eq. (|1 . 14[) was performed by 
using a Gaussian integration scheme (for details, see Ref. 6). Since each 
energy is independent of each other, the algorithm parallelizes naturally. 

It is important to find ways to benchmark this nonequilibrium DMFT 
algorithm, to ensure that it is accurate. The simplest way to do this is 
to calculate the equilibrium results within the nonequilibrium formalism 
and compare those results with the results obtained by the equilibrium 
DMFT approach. One of the most important elements is a proper choice 
of the discretization step At of the contour. These equilibrium calcula- 
tions can always help to choose the step At small enough to get cor- 
rect results (see Ref. 6). Another useful way to check the accuracy of 
the solution is to calculate the moments of the electron spectral functions 
A(t ave ,ui) — ^ k (— l/7r)ImGk(taue, <^), where t ave is the average time and 
oj is the electron frequency arising from a Fourier transform of the relative 
time (see below). We have found 7 that the lowest spectral moments in 
the Falicov-Kimball model can be calculated exactly, and they are time- 
independent even in the presence of a time-dependent electric field. In 
particular, when a spatially homogeneous time-dependent electric field is 
applied, one can find for the zeroth and first two retarded spectral moments: 




(1.18) 



/OO 
duiuA R (t ave ,cj) = -/j, + Urif = 0, 
-OO 



(1.19) 
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1 1 U 2 

dLouj 2 A R {t ave ,uj) = - + /i 2 - 2-7/zn/ + t/ 2 n/ = - + — , (1.20) 

where the second equality holds in the half-filled case. We estimate the 
accuracy of the discretization of the contour by calculating the spectral mo- 
ments and comparing them with the exact analytical results in Eqs. (|1.18|) - 
(jl.20p . 7 In general, one needs to reduce the discretization size as the inter- 
action strength increases. This is clearly seen in the equilibrium case, where 
the numerics can be well controlled because there is no dependence on the 
energy e. Surprisingly, in the presence of a field, one can use a somewhat 
larger discretization size, especially for moderate to large fields. 



1.4. Gauge invariance and physical observables 

In nonequilibrium problems, we work with two-time Green functions be- 
cause the system no longer has time-translation invariance. Wigner 33 first 
realized that it is more physical to express results in terms of average 
and relative coordinates, where the dependence on the average coordinates 
drops out in equilibrium. In our case, the relative and average times satisfy 

i ^1 ^2; tave 

while for the spatial coordinates we have 

r = ri-r 2 , r ave 

note that at this point we are restricting the time coordinates to lie on the 
real axis piece of the contour since the imaginary axis piece is not impor- 
tant for determining physical properties on the lattice (the full structure 
is only needed for the self-consistent DMFT loop, not for calculating any 
physical properties once the self-energy has been determined). We want to 
be able to convert the relative time and space coordinates into frequency 
and momentum via a Fourier transformation. Since we are working with 
a uniform electric field, we expect that the system will have no average 
spatial coordinate dependence, because it is spatially homogeneous. The 
easiest way to construct the right transformation is to create a Fourier 
transformation that makes the gauge- invariance of the problem manifest; 
the result is then called the gauge-invariant Green function, which depends 
only on the fields, not on the scalar or vector potentials. 5 The procedure 
is somewhat technical, but completely straightforward. The starting point 



h +t 2 



Ti + r 2 



(1.21) 



(1.22) 
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is a generalized Fourier transformation 

G{k,uj,r ave ,t ave ) = J d d r J dtexp[iW(k,Lu,r,t,r ave ,t ave )} 

x G\T,t,Y ave ,t ave ), (1.23) 

with W being a complicated function of its variables, in general. In equi- 
librium, when there is no external space- and time-dependent electric field, 
the Green function doesn't depend on the average coordinates t ave and 
r ave , and the transform (|1.23p is the well-known Fourier transformation 
with W(k, u>, r, t, r ave ,t ave ) = tui-r-k. 

The situation is more complicated when an external field is present. In 
this case, the field is introduced by using a specific gauge for the scalar and 
vector potential (we work with the Hamiltonian gauge) . It is important to 
have a Green function on the left hand side of Eq. (|1.23p . which doesn't 
depend on the choice of gauge so all results are manifestly independent of 
the scalar and vector potentials. Therefore, we need to construct a func- 
tion W(k,cu,r,t,r ave ,t ave ) in Eq. (|1.23p . which makes G(k,uj,r ave ,t ave ) 
invariant under the gauge transformation: 

(p (r 1 ,t 1 )^i P (v 1 ,t 1 )- dx{ 2 ,tl \ (1-24) 

A(r 1( ti) -» A(r x ,ti) + Vx(ri,*i), (1.25) 

where x( r i>£i) i s an arbitrary function. The \ function must also be used 
in the local unitary gauge transformation of the fermion operators: 

c(ri,ti) -> exp[iex(ri,ii)]c(ri,ii), (1.26) 
c + (r 2) t 2 ) -> exp[-iex(r 2 ,i 2 )]ct(r 2 ,t 2 ), (1.27) 

since it corresponds to the phase picked up by the fermions as a result of 
the local gauge transformation. Obviously, the Green function on the right 
hand side of Eq. (|1.23|) is not generically invariant in this case: 

G(ri,ii;r a ,t a ) -> exp[ie(x(ri,ti) - x(r 2 ,f 2 ))]G(ri,ti;r 2 ,t 2 ). (1-28) 

However, it is possible to show that its transform in Eq. (|1.23[) is invariant, 
when one chooses 5 

,1/2 

W(k,Lu,r,t, r ave ,t ave ) = / d\{t[u + &p(r ave + Ar, t ave + At)) 

J -1/2 

+ Ar,W + At)]}(1.29) 

(for details, see Ref. 34). 
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In the case of a spatially homogeneous electric field in the Hamiltonian 
gauge with ip(r, t) = 0, which we study in this paper, this transformation 
is 



because the function W just involves a shift of the momentum; note that 
the Green function is actually independent of r ave in this case. Hence, the 
gauge invariant Green function in the momentum representation contains 
a shift of the momentum, which depends on both the relative and average 
time coordinates. We consider the case when a constant electric field is 
turned on at time t = 0: A(t) = —~Et6(t). Then the momentum shift is 



Note that this shift does not depend on the relative time coordinate t for 
long times, t ave > \t/2\. However, in general, one has to first shift the 
momentum, and then Fourier transform the relative time to a frequency. It 
is important that the time-dependent momentum shift takes place for some 
negative average times (if the absolute value of the relative time is large 
enough, then either ti or t 2 is larger than and hence "sees" the field). 

The shift of the momentum becomes particularly simple for equal time 
Green functions, such as those needed to calculate the current flowing or 
to determine the distribution of the electrons amongst the quantum states. 
In this case, t — 0, and the momentum is shifted by — eEt ave if t ave > 
0. Therefore, gauge invariant Green functions can be obtained from the 
Hamiltonian gauge Green functions by simply shifting the momentum by 
— e~Et ave . Note that local quantities, like the local density of states or the 
local distribution function are always gauge invariant, because they are 
summed over momentum, and if the shift is the same for each momentum 
value, then we still sum over all the momentum points in the Brillouin zone. 
In cases where the relative time is nonzero, the transformation from the 
Green function in a particular gauge to the gauge-invariant Green function 
must be handled with care. Finally, one should note that in the steady state, 
where t ave — > oo, the momentum shift is also simple (— eE,t ave ); it turns out 




k ► k 6E t a ve@(tave 




(1.31) 
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that the retarded and advanced Green functions depend only on the relative 
time, but the lesser, greater, and Keldysh Green functions depend on both 
the average and relative time because there is an average-time-dependent 
shift of the momentum in Fermi-Dirac distribution functions. Caution must 
be used in trying to directly find the steady-state Green functions, because 
the Dyson equation is modified, since the momentum shift does not remove 
all average time dependence in internal variables that are integrated over 
in the GqSG term. 



1.5. Bloch electrons in infinite dimensions 

The work presented in this section is based on Ref. 31 where the original 
solution for Bloch electrons in a field was given. There the work focused 
on the Hamiltonian gauge, here we discuss the gauge-invariant formalism. 

Bloch 35 and Zener 36 originally showed that when electrons are placed 
on a perfect lattice, with no scattering, the current oscillates due to Bragg 
reflection of the wavevector as it evolves to the Brillouin-zone boundary. 
Here we show how to analyze this problem on the infinite-dimensional hy- 
percubic lattice. The noninteracting problem can be solved exactly in the 
case of an arbitrary time-dependent electric field. In particular, the nonin- 
teracting contour-ordered Green function is (in the Hamiltonian gauge 31 ): 



where /(e(k) — /i)] = l/{l+cxp[/3(e(k)— /i)]} is the Fermi-Dirac distribution 
(half-filling corresponds to \x = 0). The symbol # c (ii,i 2 ) is equal to one if 
t\ lies after ti on the contour, and is zero otherwise. Note that the Green 
function in Eq. (|1.32p is also used in the system of equations (| 1 . 14|) - (| 1 . 1 7[) 
to solve the interacting problem. 

When we have a constant electric field directed along the diagonal and 
turned on at t — 0, each component of the vector potential satisfies A(t) = 



G c k °(t!,t 2 ) 



i[/(e(k) - /x) - 8c{ti,h)] exp[i/i(ti - h)} 




(1.32) 
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—Et9(t). Then the integral that appears in the exponent of Eq. Q1.32p is 

6{-\t/2\-t ave )e{k)t (1.33) 
+ 9(-t/2-\t ave \) 

' e(k)(sin eE(t ave - t/2) + t ave + t/2) - e(k)(cos eE(t ave - t/2) - 1)" 

eE 

+ 6{t/2-\t ave \) 

' e(k)(sm eE(t ave + t/2) - t ave + t/2) + e (k) (cos eE{t ave + t/2) - 1)' 

eE 

+ 0{t ave - |t/2|)[e(k)(sine£'(W + V 2 ) 

- smeE(t ave - t/2)) + e(k)(cos eE{t ave + t/2) - coseE(t ave - t/2)) j eE 

when expressed in terms of the Wigner coordinates. 

To get the gauge-invariant Green function, we now shift the momentum 
as shown in Eq. (|1 .31|) : note that the shift is done both for the momentum 
in the exponent, and for the momentum in the Fermi-Dirac distribution. 
Two of the four cases for the exponent in Eq. (|1.33|) are easy to work out for 
the gauge-invariant Green functions. The first is the 6(— \t/2\ — t aV e) term 
which remains unchanged and the second is the 9(—t/2— \t ave \) term, which 
becomes 2e(k) sm(eEt/2). Note that both these exponents are independent 
of the average time. The average time enters for the other two terms, and 
in the argument of the Fermi-Dirac distribution. 

The retarded G^°(t 1} t2) and lesser G^°(ii,t2) Green functions can be 
obtained from Eq. (|1.32|) . by replacing the prefactor [/(e(k) — /x) — # c (£i,i 2 )] 
by —i8(ti — i 2 ) and /(e(k) — /i), correspondingly. One needs to shift the 
momentum accordingly to get the gauge-invariant retarded and lesser Green 
functions. Note that at long times the gauge-invariant retarded Green 
function depends only on relative time. 

Since the electrical current is found from the time derivative of the 
polarization operator, the current operator is determined by taking the 
commutator of the Hamiltonian (in a particular gauge) with the polariza- 
tion operator. The result, for the ath component of the current-density 
operator is 



3a(tave) — ^ ^ ^ 01 Cj c (t ave )ck(t ^ e ), (1.34) 



de(k - eA(t ave )) t 



dk a 



where we have emphasized that the operator is evaluated with a vanish- 
ing relative time (t = 0). We want the expectation value of the current 
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operator, which is found by taking expectation value of the expression in 
Eq. p.34| and noting that each component gives the same result for a 
field pointing along the diagonal. The expectation value of the product 
c\ i (t a ve)c] S i(t ave ) can be replaced by the lesser Green function G^(t ave ,0). 
So we have 



j(tave) = «e^e(k- eA(t oue ))G£(f ot , e ,0) = ie^e(k) G£(t ave ,0), 

k k 

(1.35) 



where the second equality comes from the transformation to the gauge- 
invariant Green function. The summation over momentum can be con- 
verted to a double integral over the joint density of states in Eq. (| 1 . 1 3|) . 
Substituting in the expression for the lesser Green function, yields 



The current is a periodic function of time, even though the field is time- 
independent; this effect is called a Bloch oscillation. 37 The period of the 
oscillation is equal to 2-K/eE. In order to see this oscillation in real solids, 
one needs to prepare a system where the scattering time is longer than 
the period of the Bloch oscillations. In solids, the scattering time is much 
shorter than the oscillation period, so this effect is not observed. However, 
Bloch oscillations are seen in semiconductor superlattices, where the period 
of oscillations is much shorter due to the larger lattice spacing. As we show 
in the following Sections, the effects of strong electron-electron correlations 
modify the Bloch oscillations significantly, but the driven oscillations in 
large fields survive for a surprisingly long time. 

Now we discuss the time-dependence of the density of states (DOS) for 
noninteracting electrons in a constant electric field. The DOS is found by 
using the Wigner time coordinates in Eq. (|1.21|L and making a Fourier 
transformation of the corresponding Green functions with respect to the 
relative time coordinate. In particular, the local DOS is 



Since this is a local quantity, summed over all momenta, it is automatically 
gauge-invariant. The local retarded Green function can be obtained from 




(1.36) 



where the single-particle density of states is 




(1.37) 




(1.38) 
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Eq. p.32[) . It is possible to show 31 that the steady state (t ave — > oo), for 
the case of a constant electric field turned on at t = 0, has a retarded Green 
function which satisfies 

'locitave ► OO, t) = -i6(t) ex P 



r<R 



{cos(eM) - 1} 



(1.39) 



Substitution of this expression into Eq. (|1.38p , and evaluating the Fourier 
transform with respect to the relative time, yields the steady state DOS, 
which consists of a set of delta-functions with different amplitudes (called 
the Wannier-Stark ladder 38 ). The distance between the delta-function 
peaks is equal to eE. The weight of these peaks is: 31 

r 2ir £*2 



WN 



l E 2 



ducos(Nu) exp( 



1]), 



(1.40) 



■2e 2 E 21 

for the Nth Bloch frequency, lon = eEN. It takes an infinite amount of 
time for the delta functions to develop. 
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Fig. 1.2. Density of states for noninteracting electrons with eE = 1 at different values 
of time t ave (be aware that the vertical scale changes from plot to plot). Note how the 
build up of the delta function at the Bloch frequencies is slow. 



In Fig. 11.21 we show how the DOS evolves from the time the field is 
turned on, at t ave = 0, to a large time. The DOS remains Gaussian for 
tave < — 2 and then develops large oscillations as t aV e increases. Though 
the DOS oscillates and acquires negative values in the transient regime, it 
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is possible to show (numerically) that its first three moments always satisfy 
the relevant sum rules. 



1.6. Exact solution 



In this Section we present the results for the interacting case, 6,8,9 where 
we vary the Coulomb repulsion through the metal-insulator transition that 
occurs at U c — The problem is solved by numerically solving the 

DMFT loop in Eqs. ([TT4>([l~T7|l . 

Once the Green functions and self-energy have been found by self- 
consistently solving the DMFT equations, we can extract the momentum- 
dependent lesser Green function and use Eq. (|1.35[) to find the current. 
In these calculations, we used the Green functions in a particular gauge, 
but one could easily shift to the gauge-invariant Green functions if desired. 
When the field is small, and the correlations are small, we see a damping 
of the Bloch oscillations, as expected. This is shown in the left panel of 
Fig. 11.31 One can see that the Bloch oscillations maintain their periodicity, 
but are damped as the scattering increases. As we start to approach the 
metal- insulator transition at U « 1.414, one can see the character of the 
oscillations start to change. As we move into the insulating phase, as shown 
in the right panel, the character of the oscillations changes completely, and 
we no longer see the regular Bloch structure. The oscillations seem to 
survive to much longer times than would be expected from a Boltzmann 
equation type of analysis. It remains unclear whether the steady state has 
some residual oscillations, or it goes to a constant value as predicted by 
semiclassical ideas. 




Fig. 1.3. Electric current for different values of U with fi = 10: (a) metals (U = 0, 
U = 0.25, U = 0.5, and U = 1) and (b) insulators (U = 1.5 and U = 2). 



Even more surprising is the fact that when the field is large, the current 
displays two anomalous features: (i) first, its decay is much slower than ex- 



February 6, 2008 22:40 World Scientific Review Volume - 9in x 6in paper 



18 V. Turkowski and J.K. Freericks 

pected from a semiclassical approach, where the relaxation time is inversely 
proportional to the imaginary part of the self-energy at the chemical po- 
tential, which is proportional to l/U 2 and (ii) the current develops beats 
with a beat frequency proportional to l/U. An example of this behavior 
is shown in Fig. 11.41 These beats are always present in the metallic phases 
(for large U), but disappear once one moves into the insulator. 



c 

<D 
1_ 
D 

O 



20 40 60 80 100 
Time 

Fig. 1.4. Time-dependence of the current for U = 0.5, E = 2.0, and /3 = 10. Note how 
the current has beats in its time dependence and that the decay of the current is rather 
slow. 

The time-dependence of the density of states can be calculated from 
Eq. (|1.38p . Here, we present some results for the case when the system 
is initially in the metallic phase, Fig. 11.51 What we find is that for small 
fields, the delta function peaks of the Wannier-Stark ladder get broadened, 
but the structure is still readily apparent. But as we increase the field 
strength, the behavior qualitatively changes, and in the long-time limit, 
the system evolves into a peaked structure, where the peaks are maximal 
near the edges of minibands, which are spaced apart in size by U , and the 
DOS has a local minimum in the center, where the Wannier-Stark peak 
used to appear. This is also behavior that is quite surprising. 

As the scattering increases, the DOS approaches the steady state value 
relatively quickly. This illustrates the dichotomy between the average time, 
which is important for determining the current, and the relative time, which 
is important for determining the DOS. The decay is rapid as a function of 
relative time, but is much slower as a function of average time. 

1.7. Perturbation theory 

A perturbative analysis can be performed directly on the lattice. 40 In this 
case, we do not need any DMFT loop, and we can restrict the contour to 
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Fig. 1.5. Density of states for different average times from t ave = to t aV e = 40 for 
U = 0.5, E = 1 and fi = 10. Note how the DOS develops split peaks, separated by 
U = 0.5 around the Bloch frequencies (integers here) 

be solely on the real axis. As described above, the perturbation theory is 
similar for the equilibrium and nonequilibrium cases, with the only signif- 
icant changes being that one needs to calculate with time-ordered objects 
along the contour and one needs to use noninteracting Green functions in 
the field. A strictly truncated expansion for the self-energy to second order 
in U is equal to the usual Hartree-Fock term (which vanishes at half filling) 
plus a second-order term which satisfies 

Z c W(h,t 2 ) = U 2 w 1 (l-w 1 )Gl c (h,t 2 ); (1.41) 

one can determine the retarded and lesser self-energies from this in a 
straightforward fashion, ft turns out that this truncated perturbation the- 
ory is most accurate at short times — in essence, the perturbation scries 
expansion is an expansion in a power series in time away from the time 
the field was turned on. In more conventional perturbation series in terms 
of frequency-dependent Green functions, the perturbation series is most 
accurate at high frequencies, and least accurate at low frequencies. Af- 
ter performing a Fourier transform, one can immediately see that this is 
equivalent to the perturbation theory being most accurate at short times 
and breaking down at long times. Indeed, we find that this perturbative 
treatment cannot reproduce the steady-state behavior at long times. 
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As a benchmark of our calculation, we compare the equilibrium self- 
energy found from a numerical solution of the DMFT equations to the 
perturbative result for small U. We find quite good agreement in the small- 
U range, and the equilibrium case appears to be fairly accurate up to U ~ 
0.5. For larger U values the perturbation theory breaks down — it is not 
capable of properly describing the Mott-insulating phase. 

Next we analyze the time-dependence of the electric current calculated 
by second order perturbation theory in the case when a constant electric 
field is turned on at time t = 0. Before presenting the second-order pertur- 
bative solution for the current, we briefly review the corresponding results 
from a semiclassical Boltzmann equation approximation. As was mentioned 
above, these results are qualitatively different from the exact solution. 

In the Boltzmann equation approach, one introduces a nonequilibrium 
quasiparticle distribution function / non (k, t) = —iG^(t,t), which satisfies 
the following phcnomenological equation: 

J f t ' ^ + eE(t) • V k / non (k,i) = --[/»-(k,t) - /(k)], (1.42) 
with the boundary condition: 

r°»(k,t = 0) = /(k)= — - i (1.43) 

exp[/3(e(k) - n)\ + \ 

This equation can be solved exactly (see, for example Ref. 40). Substi- 
tution of the expression for the distribution function instead of —iG < into 
Eq. (|1.35[) allows one to calculate the semiclassical current. This semiclassi- 
cal current approaches a steady state as time goes to infinity. In particular, 
in the infinite-dimensional limit one obtains: 

6 gEt f 

](t) = -V3l + eW J d ^ ef ^ 

1 - {cos{eEt) - eETsin(eEt)) e' t/r . (1.44) 

Therefore, the current is a strongly oscillating function of time for t <C t, 
and it approaches the steady-state value 



eEr 

r. 3 ^ au J — 

J u 

where 



j = --!Ljdep(e)ef(e), (1.46) 

as t/r — * oo. The steady-state current amplitude is proportional to E in 
the case of a weak field (the linear- response regime), and then becomes 
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proportional to l/E at eEr — > oo. The amplitude of the current goes to 
zero in this nonlinear regime with the field amplitude increasing. One would 
naively expect that the second-order perturbation theory would give similar 
results in the case of a weak Coulomb repulsion, since one can extract an 
effective scattering time for the equilibrium limit of the Falicov-Kimball 
model with small U: r = 1/(tt 2 U 2 ). 39 However, it will be shown below, 
that the behavior of the current calculated in second-order perturbation 
theory is rather different from the Boltzmann equation case, and closer to 
the exact numerical result at short times. 

The electric current in the second-order perturbation theory can be 
calculated by substituting the expression for the second-order lesser Green 
function into Eq. Q1.35p . In this case 



j(t) = -= de dep 2 (e,e) [ecos (eA a (t)) - esm(eA a (t))] Gf iS (t,t). 



It is difficult to find exact analytical expressions for the current, except 
for some limiting cases. Of course, in the limit U = we recover the free 
electron case result: 



where the amplitude jo of the Bloch oscillations is given by Eq. Q1.46p . 
Therefore, the general expression for the time-dependence of the electric 
current in a strictly truncated second-order perturbation theory expansion 
can be written as: 



The electric current is a superposition of an oscillating part and some other 
piece proportional to U 2 . Obviously this cannot produce a constant steady- 
state current for all small U, because the function j% is independent of U. 
This is a clear indication that the perturbation theory will hold only for 
short times. 

Numerical results for the time-dependence of the electric current calcu- 
lated from Eq. (| 1 .47)) at eE — 1 and different values of U are presented in 
Fig. 11.61 (dashed lines). Note how the current oscillates for all times within 
our finite time window. We also show the corresponding Boltzmann equa- 
tion solution and the exact solution. We use two different values for the 
Boltzmann equation — one fixes the relaxation time to the prediction from 
the equilibrium solution, while the other adjusts the relaxation time to ob- 
tain the best fit. Comparison of the perturbation theory and the Boltzmann 




(1.47) 



j(t) = j sin (eEt) , 



(1.48) 



j{t) = j sin (eEt) + U 2 j 2 {t). 



(1.49) 
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Fig. 1.6. Perturbative expansion for the electric current as a function of time for E = 
1.0, (3 = 10 and different values of U (dashed lines). The solid and dotted lines correspond 
to the exact DMFT solution and the Boltzmann equation (BE) solution, respectively. 
The dash-dotted lines (BE2) are the Boltzmann equation result with a phenomenological 
relaxation time r = a/ '(ir 2 U 2 ), while a = 20 in Figs, a) and b), and a = 36.5 in Fig. c). 



equation solution shows that they are close at short times, but at longer 
times the PT current remains oscillating, while the Boltzmann equation 
solution approaches a steady state. Moreover, at times longer than ~ 2/U 
the perturbation theory breaks down showing an oscillating current with 
increasing amplitude. At times shorter than 2/U the perturbation theory 
solution is close to the exact result, displaying an oscillating current with 
decreasing amplitude. It is also possible to fit the Boltzmann equation 
results to the exact and PT solution at short times if one chooses the re- 
laxation time t — a/ (tt 2 U 2 ), where a ~ 20 — 30, which is much larger than 
a = 1 in the case of the second order perturbation theory. 39 These results 
clearly show that the semiclassical approach, with one effective time vari- 
able that is damped on the timescale of the relaxation time is not sufficient 
to describe the behavior in the quantum case. 

The perturbation theory calculations at different values of the electric 
field give results similar to the results presented in Fig. 11.61 Numerical 
analysis shows that the agreement between the perturbation theory results 
and the exact results is better when eE is larger than U. In fact, it is 
possible to find an analytical expression for the current in the case of large 
electric fields: 



j{t) ~ -— / dep(e)ef(e) 



1 - U 2 B(f3) - ^-t 2 



sm(eEt), (1.50) 



where B{(3) is a positive decreasing function of temperature: 40 0.25 < 
B(J3) < 0.5. 



February 6, 2008 22:40 World Scientific Review Volume - 9in x 6in paper 



Nonequilibrium dynamical mean-field theory of strongly correlated electrons 23 

1.8. Conclusions 

To conclude, we have presented some results on the nonequilibrium prop- 
erties of the Falicov-Kimball model of strongly correlated electrons in the 
limit of infinite dimensions. Despite the simplicity of the model, the so- 
lutions show that strong electron-electron correlations result in nontrivial 
behavior. The dynamical mean-field theory approximation is believed to be 
a precise method to solve strongly correlated problems in three dimensions. 
Therefore, we believe that some of the results, like the long-time Bloch 
oscillations of the current, beats in the current for strong fields and the 
splitting of the Wanier-Stark peaks could be observed in bulk systems with 
dominant electron-electron scattering in the presence of a strong electric 
field. Such a field can be present, in particular, in nanostructures, where 
a moderate external electric potential can produce strong (uniform) elec- 
tric fields due to the small size of the systems. One might also be able 
to observe this behavior in mixtures of heavy and light atoms trapped in 
optical lattices. In addition, we demonstrated that the perturbation theory 
solution cannot be used to study the long-time behavior of the system. It 
would be interesting to generalize these results to more complicated models 
and to lower dimensions, where we expect qualitatively similar behavior. 
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